#R Replication file for Macropartisanship Revisited
#Green, Hamel, and Miller, Perspectives 2023

#Prior to running this script, load "MPID R Replication.RData"


#Packages
library(tidyverse)
library(strucchange)
library(TSstudio)
library(forecast)
library(ggseas)
library(plotly)
library(stargazer)
library(lubridate)
library(cowplot)
library(texreg)

#Replication of Paper Results

#Table 1: Descriptive Statistics
descriptives<-macro.tibble %>%
  select(macropid_adj,party_ics,party_approval,ics_unadjusted,gallup_approval) %>%
  as.data.frame()
summary(descriptives)

#Figure 1
ggplot(macro.tibble, aes(x=Time, y=macropid_adj)) +
  geom_line() + 
  geom_vline(xintercept = 1988, linetype="dashed", 
             color = "red", size=1) +
  xlab("") +
  ylab("Macropartisanship") +
  theme_bw() +
  annotate(geom="text", x=1970, y=.55, label="MES Period",
           color="black") +
  annotate(geom="text", x=2005, y=.65, label="Modern Period",
           color="black") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))

#Figure 2
long<-macro.tibble %>%
  select(Time, macropid_adj, ics, gallup_approval) %>%
  mutate(ics=ics/100) %>%
  rename('Consumer Sentiment' = ics,
         'Macropartisanship' = macropid_adj,
         'Presidential Approval' = gallup_approval) %>%
  pivot_longer(!Time, names_to = "indicator", values_to =  "value")

ggplot(long, aes(x = Time, y = value, color = indicator, linetype = indicator)) +
  geom_line(size=.8) +
  scale_color_manual(values=c('Green','Blue','Red')) +
  scale_linetype_manual(values=c("solid", "dotdash", "dotted")) +
  geom_vline(xintercept = 1961, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1963, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1969, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1974.3, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1977, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1981, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1989, 
             color = "black", size=.5) +
  geom_vline(xintercept = 1993, 
             color = "black", size=.5) +
  geom_vline(xintercept = 2001, 
             color = "black", size=.5) +
  geom_vline(xintercept = 2009, 
             color = "black", size=.5) +
  geom_vline(xintercept = 2017, 
             color = "black", size=.5) +
  geom_vline(xintercept = 2021, 
             color = "black", size=.5) +
  xlab("") +
  ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        plot.margin = margin(.5, .55, 0, 0, "cm"),
        legend.position="bottom",
        legend.title=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x = element_blank()) +
  annotate(geom="text", x=1957, y=.15, label="Eisenhower",
           color="black", size=3.25) +
  annotate(geom="text", x=1962, y=.15, label="JFK",
           color="black", size=3.25) +
  annotate(geom="text", x=1966, y=.15, label="LBJ",
           color="black", size=3.25) +
  annotate(geom="text", x=1971.5, y=.15, label="Nixon",
           color="black", size=3.25) +
  annotate(geom="text", x=1975.5, y=.15, label="Ford",
           color="black", size=3.25) +
  annotate(geom="text", x=1979, y=.15, label="Carter",
           color="black", size=3.25) +
  annotate(geom="text", x=1985, y=.15, label="Reagan",
           color="black", size=3.25) +
  annotate(geom="text", x=1991, y=.15, label="HW Bush",
           color="black", size=3.25) +
  annotate(geom="text", x=1997, y=.15, label="Clinton",
           color="black", size=3.25) +
  annotate(geom="text", x=2005, y=.15, label="W Bush",
           color="black", size=3.25) +
  annotate(geom="text", x=2013, y=.15, label="Obama",
           color="black", size=3.25) +
  annotate(geom="text", x=2019, y=.15, label="Trump",
           color="black", size=3.25) +
  coord_cartesian(ylim = c(.25,1.05), clip = "off") +
  scale_x_continuous(limits=c(1953,2021), expand=c(0,0), breaks = scales::pretty_breaks(n = 10))


##### MODELS

#Table 2

table.2.data<-macro.tibble %>%
  mutate(lag.mpid=lag(macropid_adj, k = 1),
         modern=year>1987) 

#ICS Only

#Fit Models for all three periods

#Table 2, Col 1
col.1<-lm(macropid_adj~lag.mpid+party_ics+party_ind, data=table.2.data)

#Table 2, Col 2
col.2<-lm(macropid_adj~lag.mpid+party_ics+party_ind, data=subset(table.2.data, modern==FALSE))

#Table 2, Col 3
col.3<-lm(macropid_adj~lag.mpid+party_ics+party_ind, data=subset(table.2.data, modern==TRUE))


#Approval Only

#Table 2, Col 4
col.4<-lm(macropid_adj~lag.mpid+party_approval+party_ind, data=table.2.data)

#Table 2, Col 5
col.5<-lm(macropid_adj~lag.mpid+party_approval+party_ind, data=subset(table.2.data, modern==FALSE))

#Table 2, Col 6
col.6<-lm(macropid_adj~lag.mpid+party_approval+party_ind, data=subset(table.2.data, modern==TRUE))


#Compile Table 2 in .tex
texreg(list(col.1, col.2, col.3, col.4, col.5, col.6),
        digits=4)



####Breakpoint Analysis: Locating breaks in the "sc" column of Table 3

macro.model.1 <- macropid_adj ~ party_ics + party_ind + lag_macro
macro.model.2 <- macropid_adj ~ party_approval + party_ind + lag_macro


#Find breakpoints:
break.1 <- breakpoints(macro.model.1, 
                       data = macro.pid, 
                       h = .2,
                       breaks=4)
#Specification 1 Breakpoints, allowing for up to 4:
summary(break.1)


break.2 <- breakpoints(macro.model.2, 
                       data = macro.pid, 
                       h = .2,
                       breaks=4)
#Specification 2 Breakpoints, allowing for 4:
summary(break.2)

macro.tibble<-macro.tibble %>%
  mutate(lag_approval=lag(party_approval, k=1))

#Appendix Tables:
#1 and 2 done, need 3-5


#####Table A1

#(1,0,1) Models
#ICS and Approval

actual<-macro.tibble %>%
  select(macropid_adj) %>%
  ts(start=c(1953, 1), end=c(2021, 1), frequency=4)

full.covs<-macro.tibble %>%
  select(party_ind,party_approval, party_ics) %>%
  ts(start=c(1953, 1), end=c(2021, 1), frequency=4)

train<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(1:140) %>%
  ts(start=c(1953, 1), end=c(1987, 4), frequency=4)

train.covs<-macro.tibble %>%
  select(party_ind,party_approval, party_ics) %>%
  slice(1:140) %>%
  ts(start=c(1953, 1), end=c(1987, 4), frequency=4)

test<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(141:273) %>%
  ts(start=c(1988, 1), end=c(2021, 1), frequency=4)

test.covs<-macro.tibble %>%
  select(party_ind,party_approval, party_ics) %>%
  slice(141:273) %>%
  ts(start=c(1988, 1), end=c(2021, 1), frequency=4)

#Fit Models for all three periods

#Table A1, Column 1
armaapp.full<-Arima(actual, order=c(1,0,1), xreg=full.covs)

#Table A1, Column 2
armaapp.train<-Arima(train, order=c(1,0,1), xreg=train.covs)

#Table A1, Column 3
armaapp.test<-Arima(test, order=c(1,0,1), xreg=test.covs)


#(0,0,1) Lagged DV Models

full.covs<-macro.tibble %>%
  select(party_ind,party_approval, party_ics)
mp<-macro.tibble$macropid_adj
full.covs$lag.mpid<-lag(mp, k = 1)
full<-full.covs %>%
  slice(2:273) 

full.covs<-full %>%
  ts(start=c(1953, 2), end=c(2021, 1), frequency=4)

actual<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(2:273) %>%
  ts(start=c(1953, 2), end=c(2021, 1), frequency=4)

train.covs<-full %>%
  slice(1:139) %>%
  ts(start=c(1953, 2), end=c(1987, 4), frequency=4)

train<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(2:140) %>%
  ts(start=c(1953, 2), end=c(1987, 4), frequency=4)

test.covs<-full %>%
  slice(140:272) %>%
  ts(start=c(1988,1), end=c(2021,1), frequency=4)

test<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(141:273) %>%
  ts(start=c(1988, 1), end=c(2021, 1), frequency=4)

#Fit Models for all three periods

#Table A1, Column 4
armaapplag.full<-Arima(actual, order=c(0,0,1), xreg=full.covs)

#Table A1, Column 5
armaapplag.train<-Arima(train, order=c(0,0,1), xreg=train.covs)

#Table A1, Column 6
armaapplag.test<-Arima(test, order=c(0,0,1), xreg=test.covs)


#####Table A2

#MA(1) Models With Lagged DV

#ICS Only

full.covs<-macro.tibble %>%
  select(party_ics,party_ind)
mp<-macro.tibble$macropid_adj
full.covs$lag.mpid<-lag(mp, k = 1)
full<-full.covs %>%
  slice(2:273) 

full.covs<-full %>%
  ts(start=c(1953, 2), end=c(2021, 1), frequency=4)

actual<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(2:273) %>%
  ts(start=c(1953, 2), end=c(2021, 1), frequency=4)

train.covs<-full %>%
  slice(1:139) %>%
  ts(start=c(1953, 2), end=c(1987, 4), frequency=4)

train<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(2:140) %>%
  ts(start=c(1953, 2), end=c(1987, 4), frequency=4)

test.covs<-full %>%
  slice(140:272) %>%
  ts(start=c(1988,1), end=c(2021,1), frequency=4)

test<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(141:273) %>%
  ts(start=c(1988, 1), end=c(2021, 1), frequency=4)

#Fit Models for all three periods

#Table A2, Col 1
arma4.full<-Arima(actual, order=c(0,0,1), xreg=full.covs)

#Table A2, Col 2
arma4.train<-Arima(train, order=c(0,0,1), xreg=train.covs)

#Table A2, Col 3
arma4.test<-Arima(test, order=c(0,0,1), xreg=test.covs)


#Approval Only

full.covs<-macro.tibble %>%
  select(party_approval, party_ind)
mp<-macro.tibble$macropid_adj
full.covs$lag.mpid<-lag(mp, k = 1)
full<-full.covs %>%
  slice(2:273) 

full.covs<-full %>%
  ts(start=c(1953, 2), end=c(2021, 1), frequency=4)

actual<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(2:273) %>%
  ts(start=c(1953, 2), end=c(2021, 1), frequency=4)

train.covs<-full %>%
  slice(1:139) %>%
  ts(start=c(1953, 2), end=c(1987, 4), frequency=4)

train<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(2:140) %>%
  ts(start=c(1953, 2), end=c(1987, 4), frequency=4)

test.covs<-full %>%
  slice(140:272) %>%
  ts(start=c(1988,1), end=c(2021,1), frequency=4)

test<-macro.tibble %>%
  select(macropid_adj) %>%
  slice(141:273) %>%
  ts(start=c(1988, 1), end=c(2021, 1), frequency=4)

#Fit Models for all three periods

#Table A2, Col 4
arma5.full<-Arima(actual, order=c(0,0,1), xreg=full.covs)

#Table A2, Col 5
arma5.train<-Arima(train, order=c(0,0,1), xreg=train.covs)

#Table A2, Col 6
arma5.test<-Arima(test, order=c(0,0,1), xreg=test.covs)

#Table A3:

#ICS Only

#Table A3, Col 1
a3.col.1<-lm(macropid_adj~lag.mpid+party_ics+party_ind +
            kennedy1 + johnson1 + nixon1 +
          ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
          trump1 + biden1, data=table.2.data)

#Table A3, Col 2
a3.col.2<-lm(macropid_adj~lag.mpid+party_ics+party_ind +
            kennedy1 + johnson1 + nixon1 +
            ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
            trump1 + biden1, data=subset(table.2.data, modern==FALSE))

#Table A3, Col 3
a3.col.3<-lm(macropid_adj~lag.mpid+party_ics+party_ind +
            kennedy1 + johnson1 + nixon1 +
            ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
            trump1 + biden1, data=subset(table.2.data, modern==TRUE))


#Approval Only

#Table A3, Col 4
a3.col.4<-lm(macropid_adj~lag.mpid+party_approval+party_ind +
            kennedy1 + johnson1 + nixon1 +
            ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
            trump1 + biden1, data=table.2.data)

#Table A3, Col 5
a3.col.5<-lm(macropid_adj~lag.mpid+party_approval+party_ind +
            kennedy1 + johnson1 + nixon1 +
            ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
            trump1 + biden1, data=subset(table.2.data, modern==FALSE))

#Table A3, Col 6
a3.col.6<-lm(macropid_adj~lag.mpid+party_approval+party_ind +
            kennedy1 + johnson1 + nixon1 +
            ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
            trump1 + biden1, data=subset(table.2.data, modern==TRUE))


#Compile Table A3 in .tex
texreg(list(a3.col.1, a3.col.2, a3.col.3, a3.col.4, a3.col.5, a3.col.6),
       digits=4)


#Table A4

#Table A4, Col 1
a4.col.1<-lm(macropid_adj~lag.mpid+party_ics+party_approval+party_ind +
               kennedy1 + johnson1 + nixon1 +
               ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
               trump1 + biden1, data=table.2.data)

#Table A4, Col 1
a4.col.2<-lm(macropid_adj~lag.mpid+party_ics+party_approval+party_ind +
               kennedy1 + johnson1 + nixon1 +
               ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
               trump1 + biden1, data=subset(table.2.data, modern==FALSE))

#Table A4, Col 1
a4.col.3<-lm(macropid_adj~lag.mpid+party_ics+party_approval+party_ind +
               kennedy1 + johnson1 + nixon1 +
               ford1 + carter1 + reagan1 + hwbush1 + clinton1 + wbush1 + obama1 +
               trump1 + biden1, data=subset(table.2.data, modern==TRUE))


#Compile Table A4 in .tex
texreg(list(a4.col.1, a4.col.2, a4.col.3),
       digits=4)



#Appendix Table 5
#Comparing Apprvoal to Lagged Approval in the Full Series

#Reproduce the Model from Table 2, Col 4
lag.approval<-lm(macropid_adj~lag.mpid+lag_approval+party_ind, data=table.2.data)

approval.potus<-lm(macropid_adj~lag.mpid+party_approval+party_ind
                   + kennedy1 + johnson1 + nixon1 + ford1+ carter1+
                     reagan1+ hwbush1 + clinton1 +wbush1 + obama1+
                     trump1 + biden1, data=table.2.data)

lag.approval.potus<-lm(macropid_adj~lag.mpid+lag_approval+party_ind
                       + kennedy1 + johnson1 + nixon1 + ford1+ carter1+
                         reagan1+ hwbush1 + clinton1 +wbush1 + obama1+
                         trump1 + biden1, data=table.2.data)


#Compile Table in .tex
#Results show much smaller coefficient for lagged approval than approval
texreg(list(col.4, lag.approval, approval.potus, lag.approval.potus),
       digits=4)
